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ABSTRACT 

We present analysis of the evolution of dark matter halos in dense environments of groups and clusters 
in dissipationless cosmological simulations. The premature destruction of halos in such environments, 
known as the overmerging, reduces the predictive power of A-body simulations and makes difficult 
any comparison between models and observations. We analyze the possible processes that cause the 
overmerging and assess the extent to which this problem can be cured with current computer resources 
and codes. Using both analytic estimates and high resolution numerical simulations, we argue that 
the overmerging is mainly due to the lack of numerical resolution. We find that the force and mass 
resolution required for a simulated halo to survive in galaxy groups and clusters is extremely high and 
was almost never reached before: ~ 1 — 3 kpc and 10 8 — 10 9 Mq , respectively. We use the high-resolution 
Adaptive Refinement Tree (ART) A-body code to run cosmological simulations with the particle mass 
of w 2 x lO 8 /!^ 1 Mq and the spatial resolution of « 1- 2/i _1 kpc, and show that in these simulations 
the halos do survive in regions that would appear overmerged with lower force resolution. Nevertheless, 
the halo identification in very dense environments remains a challenge even with the resolution this high. 
We present two new halo finding algorithms developed to identify both isolated and satellite halos that 
are stable (existed at previous moments) and gravitationally bound. 

To illustrate the use of the satellite halos that survive the overmerging, we present a series of halo 
statistics, that can be compared with those of observed galaxies. Particularly, we find that, on average, 
halos in groups have the same velocity dispersion as the dark matter particles, i.e. do not exhibit 
significant velocity bias. The small-scale (100 kpc - 1 Mpc) halo correlation function in both models 
is well described by the power law £ oc r~ and is in good agreement with observations. It is slightly 
antibiased (6 « 0.7—0.9) relative to the dark matter. To test other galaxy statistics, we use the maximum 
of halo rotation velocity and the Tully-Fisher relation to assign the luminosity to the halos. For two 
cosmological models, a flat model with the cosmological constant and fin = 1 — ^a = 0.3, h = 0.7, and 
a model with a mixture of cold and hot dark matter and Slo = 1.0, fi„ = 0.2, h — 0.5, we construct 
the luminosity functions and evaluate mass-to-light ratios in groups. Both models produce luminosity 
functions and the mass-to-light ratios (~ 200 — 400) that are in a reasonable agreement with observations. 
The latter implies that the mass-to-light ratio in galaxy groups (at least for M V i r <J 3 x 10 13 /i -1 M Q 
analyzed here) is not a good indicator of fin. 

Subject headings: cosmology: theory - large-scale structure of universe - methods: numerical 



1. INTRODUCTION 

It is generally believed that the dark matter (DM) con- 
stitutes a large fraction of the mass in the Universe and 
thus significantly affects, and on some scales dominates, 
the process of galaxy formation. Observational evidence 
for large fractions of DM in galaxies, groups, and galaxy 
clusters ranges from flat rotational curves of spiral (e.g., 
Faber & Gallagher, 1979; Rubin et al. 1985; Persic, 
Salucci, & Stel 1996; Courteau & Rix, 1997) and X-ray 
emission and mass-to-light ratios of elliptical (Forman et 
al.1985; Rix 1997; Brighenti & Mathews 1997) galaxies to 
the baryon fractions in clusters of galaxies (White et al. 
1993; Evrard 1997). For galaxies, the extent of the DM 
halos estimated using satellite dynamics, is ~ 0.2 — 0.5ft -1 



Mpc 2 (Zaritsky & White 1994; Carignan et al.1997). A 
convincing evidence for substantial amounts of dark mat- 
ter even in the very inner regions of galaxies comes from 
the recent HI studies of the dwarf and low surface bright- 
ness (LSB) galaxies. The observed amounts of stars and 
gas in some of these galaxies can account for less than 10% 
of the observed rotational velocities at the last measured 
point of the rotation curves, (e.g., Carignan & Freeman 
1988; Martimbeau, Carignan, & Roy 1994; de Blok & Mc- 
Gaugh 1997). 

If the observed galaxies have large DM halos, then N- 
body simulations can, in principle, be used to predict dis- 
tribution of the dark matter component, to associate the 
simulated DM halos with galaxies, and to predict the bulk 
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properties of these galaxies such as position, mass, and 
size. One should be able then to make predictions about 
spatial distribution and motion of these simulated galaxies 
and compare these predictions with corresponding obser- 
vations. Unfortunately, the dissipationless numerical sim- 
ulations have been consistently failing to produce galaxy- 
size dark matter halos in dense environments typical for 
galaxy groups and clusters (e.g., White 1976; van Kampen 
1995; Summers, Davis, & Evrard 1995; Moore, Katz, & 
Lake 1996). This apparent absence of substructure in the 
virialized objects, known as the overmerging problem, re- 
flects the fact that simulated galaxies seem to merge much 
more efficiently in comparison with real galaxies in groups 
and clusters. In the central regions of a cluster (<~ 500 
kpc ), the "overmerging" erases not only large-scale sub- 
structure, but also any trace of small halos that could be 
associated with "galaxies" 3 , leaving a smooth giant lump 
of dark matter. 

The overmerging problem was traditionally explained by 
the lack of dissipation in A-body simulations (e.g., Katz, 
Hernquist & Weinberg 1992; Summers et al. 1995). In- 
deed, the DM halos are much larger than baryonic extent 
of the galaxies due to the dissipational nature of the lat- 
ter. The radiative cooling, for example, allows baryonic 
component to sink into the center of the DM halo where 
it forms a compact, tightly bound object. In dense envi- 
ronments the large DM halo can be easily stripped by the 
tidal field of a galaxy cluster or group, whereas the more 
compact and denser gas clump may survive (Summers et 
al. 1995). Although it is clear that to produce a realis- 
tic galaxy we need to include the energy dissipation by 
baryons, it is not clear whether the dissipation is vital for 
the halo survival in a cluster. 

Two arguments can be presented against the traditional 
explanation for the overmerging. First, if the dissipation 
helps galaxies to survive in clusters, then galaxies should 
be dominated by baryons at all scales within their visible 
extent. Most of the observed galaxies, however, appear 
to have a substantial fraction of DM inside their optical 
radius (e.g., Persic et al. 1996). The survival of a galaxy 
dominated at its optical radius by the dark matter will de- 
pend mostly on the dark matter, not on the baryons. The 
DM dominated dwarfs must have been tidally disrupted, 
but dwarf galaxies are observed in clusters (e.g. Smith, 
Driver, & Phillipps 1997; Lopez-Cruz et al. 1997). Sec- 
ond, as we argue in this paper, even in the absence of the 
baryons DM halos are dense enough to survive inside clus- 
ters and to be identified, provided that simulations have 
sufficient resolution in both mass and force (similar con- 
clusion was reached by Moore et al. 1996). 

The main goal of this paper is to demonstrate that 
with a sufficient computational effort the overmerging 
problem can be tamed (at least to some extent) even 
in the purely dissipationless simulations. The computa- 
tional costs are higher than cost of an average cosmological 
A-body simulation. However, they may be considerably 
lower than computational expense of the corresponding A- 
body+hydro simulations. This especially true in the case 



of large- volume (~ 50 — 100ft- -1 Mpc) simulations required 
to quantify the statistical properties of the "galactic" pop- 
ulations such as correlation function, pairwise velocity dis- 
persions etc. To understand how the problem should be 
dealt with, it is important to understand what processes 
lead to the overmerging. 

Several recent studies have addressed this question from 
different viewpoints and using different numerical and ana- 
lytical techniques. Thus, for example, van Kampen (1995) 
studied formation of galaxy clusters in purely dissipation- 
less simulations and concluded that two-body relaxation 
and tidal disruption are primarily responsible for the over- 
merging. He found, however, that the two-body effects are 
important only for the smallest halos (<; 30 DM particles), 
in quantitative agreement with experiments of Moore et al. 
(1996). The latter study addressed also effects of particle- 
halo and halo-halo heating on the survival of DM halos. 
The authors concluded that particle-halo heating does not 
pose a problem as long as DM particle mass is <; 1O 1O M0, 
but halo-halo heating may be important if force resolu- 
tion is not adequate (;> lOkpc). The general conclusion 
of Moore et al. is that the overmerging problem is due 
mainly to tidal heating by the cluster and halo-halo heat- 
ing, both effects being enhanced by poor force resolution. 
They note also that these effects depend crucially on the 
density structure of DM halos. 

The density profile of dark matter halos in the CDM 
models is now known reasonably well. Navarro, Frcnk & 
White (1995, 1996, 1997, hereafter NFW) gave an analyt- 
ical fit that describes with a reasonable accuracy the den- 
sity profiles of DM halos formed in the standard cold dark 
matter scenario over large range of scales and masses. The 
analytical form of the density profile advocated by NFW, 
allows one to estimate tidal effects and effects of dynamical 
friction analytically. We present these estimates in §2. 

Two recent numerical studies (Tormen, Diaferio & Syer 
1998; and Ghigna et al. 1998) present evidence that higher 
resolution significantly aleviates the overmerging problem 
in dissipationless simulations. Both studies conclude that 
many halos survive for a long time after falling onto a large 
halo, although overmerging problem persists to a certain 
degree in the central dense region of the large cluster-size 
halo. Tormen et al. also point out that there are real 
physical processes which would lead to the erase of sub- 
structure even in the ideal very high-resolution simulation. 
Namely, the dynamical friction drives massive satellites to 
the cluster core where they get tidally stripped and quickly 
disrupted. This effect is, of course, real and should be 
distinguished from artificial overmerging caused by insuf- 
ficient numerical resolution. We will address these issues 
in §2. 

Two goals of the study presented in this paper are 1) 
to make an approximate estimates of the effects leading 
to the overmerging for the halos with the NFW density 
distribution; and 2) to demonstrate that dissipationless 
simulations with sufficiently high resolution in force and 
mass are affected by the overmerging to a considerably 
lesser degree. In other words, we present an attempt to 
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or "galaxy-size halo" to indicate the DM halos formed in the simulations which could be associated with places where luminous baryons could 
reside. 
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estimate to which extent the overmerging problem can be 
solved with numerical resolution that can be achieved with 
current codes and computational resources. The goal is 
worthwhile. If the problem can be minimized at a com- 
putational cost that is not prohibitive, the dissipationlcss 
simulations can be used for a direct study of the statistical 
properties of "galaxies" . This may include studies of their 
spatial distribution, velocity field, environmental effects 
etc. This may allow us to make a step towards solution of 
the long-standing and particularly important issue of the 
galactic bias. 

The overall plan of the paper is as follows. In §2 we 
discuss the numerical and physical effects which may lead 
to the erasure of substructure inside the dense massive ha- 
los. Specifically, we present analytical estimates for tidal 
disruption of dark matter halos and effects of dynamical 
friction assuming that halos are described by the NFW 
density profile. We use these calculations to make a rough 
estimate of what resolution would be required to mini- 
mize the overmerging. Numerical simulations and cos- 
mological models are discussed in §3. In §4 we discuss 
difficulties associated with identification of dark matter 
halos in very dense regions. We describe two new halo 
finding algorithms developed to handle the halos in such 
environments. In §5 we present results of high-resolution 
dissipationless simulations, illustrate the performance of 
the new halo finding algorithms and discuss the degree 
to which the simulations arc affected by the overmerging. 
Our conclusions are illustrated using well-known statistics 
such as two-point correlation function, velocity bias, lu- 
minosity function and (M/L) ratio. The conclusions are 
summarized in §6. 

2. SURVIVAL OF HALOS IN CLUSTERS 

2.1. Numerical effects vs. physical effects 

As was noted in the introduction, there is a number 
of processes potentially contributing to the erasure of sub- 
structure in clusters and groups (van Kampen 1995; Moore 
et al. 1996). Some of these processes are due to numerical 
limitations of a simulation, while others are real physical 
effects. 

The major effects caused by numerical limitations are 
particle evaporation due to the two-body relaxation (e.g., 
Carlberg 1994; van Kampen 1995; Moore et al. 1996), 
particle-halo heating (Carlberg 1994; Moore et al. 1996), 
and premature tidal disruption due to an insufficient force 
resolution. The two-body evaporation in an iV-body sys- 
tem is a well- understood process (e.g., Binney & Tremaine 
1987). However, it was shown (van Kampen 1995; Moore 
et al. 1996) that the two-body evaporation is important 
only for the halos with <; 30 particles. Thus, if a sim- 
ulation has a particle mass of <; 10 9 /i _1 M Q , this effect 
is not important for halos in the mass range of interest 
(;> lO 11 /^ 1 Mq). It may become important, however, if 
halo looses most of its mass due to the tidal stripping. The 
particle-halo heating was suggested by Carlberg (1994) as 
an explanation for overmerging in his simulation. How- 
ever, Moore et al. (1996) showed that the time-scale 
for this process is large for typical numerical parameters. 
Finally, gravitational softening, r so f t , imposed in an TV- 
body simulation usually leads to a constant density core 
r c rj r so f t in the halo center. Poor force resolution may, 



thus, result in halos with artificially large cores which will 
be tidally disrupted faster than if they would have a more 
compact density distribution. 

The real physical processes which may lead to the era- 
sure of substructure include dynamical friction, tidal strip- 
ping, and halo-halo heating. The dynamical friction drives 
DM halos towards the high-density cluster center where 
they get tidally stripped and merge with the central mas- 
sive object. Dynamical friction can be important for some 
halos and is probably responsible for the presence of mas- 
sive central cD galaxies in observed clusters (e.g., Merritt 
1985). The importance of dynamical friction for each par- 
ticular halo will depend on the halo mass and details of its 
orbit. We will estimate these dependencies in §2.4. Tidal 
stripping occurs simply because at some distance from the 
center the tidal force from a cluster is stronger than grav- 
itational force of a parent halo and particles beyond this 
radius become unbound from a halo and dissolve in the 
ambient diffuse medium. This effect will be discussed in 
the next section. Finally, Moore et al. (1996) argued 
that heating due to the close passages of DM halos on 
their orbits may lead to significant mass losses. Their esti- 
mate, however, was based on experiment in which cluster 
is static. In real clusters only a few halos will be present for 
its entire lifetime, the bulk of the halos being accreted over 
period of the cluster evolution. This process is probably 
less important than the more efficient tidal stripping. 

Numerical and physical effects are, of course, closely re- 
lated. Tidal force of the cluster and close encounters of 
individual dark halos result in effective stripping of the 
peripheral parts of the halos. This is a real physical effect. 
But after the stripping is done, two-body relaxation can 
become artificially short and can result in evaporation of 
halos if the number of particles in a halo is too small. Due 
to the dynamical friction halos tend to spiral down towards 
the center of the cluster where they merge with the cen- 
tral halo. This is a real effect. But it can be exacerbated 
by artificially high energy dissipation in hydrodynamical 
simulations resulting in too compact halos. 

It is clear that physical processes will operate and will 
tend to erase substructure even in an ideal simulations 
of infinite resolution. The numerical effects, on the other 
hand, can be cured by improving the spatial and mass res- 
olution of the simulations and/or by inclusion of additional 
physics (e.g., missing gas dissipation). It is the numerical 
effects that we will focus on in this paper. Throughout the 
rest of the paper, we will be using the term "overmerging 
problem" only for processes related to the numerical ef- 
fects, which appear due to the lack of resolution or due to 
missing physical effects. 

2.2. Tidal disruption of halos: analytical estimates 

The main reason for the erasure of the substructure in 
clusters is the tidal interaction of individual halos with 
the cluster potential. This can be now estimated reli- 
ably, without assuming that halos are isothermal spheres 
or have King profiles, as was typically the case in the past 
(Moore et al. 1996). The density profiles of dark matter 
halos in the CDM models is now known reasonably well 
for a large range of masses and for a variety of cosmo- 
logical models (NFW). Below we give analytic estimates 
for various halo properties in the standard f2 = 1 CDM 
model. The low-density flat CDM model with cosmolog- 
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ical constant (ACDM; e.g., Carrol, Press & Turner 1992) 
model predicts similar profiles, if the the overdensity of a 
collapsed object is adjusted properly to take into account 
the change in ft (Lahav et al. 1991; Eke et al. 1996). 
The NFW density profile is given by 



P(r)= /fV M(r)=M vir .|M 
r(r + r s y f(C) 

/(x)=ln(l+x) 



(1) 



1 + .T 



r 



where r s and po are the characteristic radius and density 
of the halo, M v ; r is the virial mass, r v i r is the virial radius, 
and C is the concentration for a halo defined as follows: 



C 



rv ,^) = 443>,-Vc( AW1 °7" Ma )' /3 .(2) 

Here, p cr is the critical density of the Universe and (5th is 
the overdensity (<5/9/p ma ttcr) of a collapsed object accord- 
ing to the top-hat model of spherical collapse. For the 
CDM model 5 t h ~ 200. Note that our definition of r v ; r 
differs from that used by NFW, who use <5 t h = 200 for 
all cosmological models. We use values predicted by the 
top-hat model which gives, for example, <5 t h ~ 340 for the 
ACDM model with fi = 0.3. 

The concentration C is a function of mass M v - U . In the 
case of the CDM model 



C w 124(M vir /lft _1 M Q ) 



-0.084 



(3) 



Typical values for the concentration C range from C ~ 12 
for M vir = lO 12 /^ 1 M Q to C w 7 for M vir = lO 15 /^ 1 M Q 
(NFW). Using these definitions wc can write mass M(r), 
orbital frequency w(r), and gravitational potential 4>{r): 



2 GM GM vir /(or) 
w ( r ) = — =- 



<f>(r) 



rf/(C) x 3 ' 
GM vir /W + ifj 



r./(C) 



a; 



(4) 



In some cases it is more convenient to define properties 
of halos by using maximum rotational velocity V max = 
(GM/r) \ max , halo parameter which is probably most 
easily related to an observable quantity, instead of the con- 
centration C or the virial mass M vlr . For profile Eq.(l) the 
maximum of the rotational velocity occurs at r max w 2r s . 
This gives 



V = 

r max 



M(r) = 



GM V 



/(2) 
2/(C) : 



r s V max 2f(x) 



G f(2) 
\ /I = -2<f>(r) = W r 



/(2) w 0.432 
V 2 



(5) 
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where V esc is the escape velocity at the distance r from 
the cluster center. 



The tidal radius, r t , of a small halo with mass m and 
maximum rotational velocity v max moving at a radius R 
from the center of a large halo with mass M{R) and V max , 
is the minimum of two radii: (1) a radius at which the 
gravity force of the small halo F grav is equal to the tidal 
force of the large halo F t ide, and (2) a radius defined by 
the resonances between the force the small halo exerts on 
the particle and the tidal force by the large halo. If r 
is the distance of a particle from the center of the small 
halo, then the condition F grav (r) — Fud e (r; R) gives an 
equation for the tidal radius r t : 



m(r t ) 
M(R) 



= 2- 



R dM 

MM' 



f{Xy) 
/(**) 



X r \ 

xr) 



3 / v \ 2 
Rs ^ma. 



(l + x R )\f{x R ) 



(6) 

where x r = r t /r s and xr = R/R s . The last equation can 
be solved numerically. 

It was argued (e.g., Weinberg 1994ab, 1997) that ef- 
fective tidal stripping can occur at smaller radius defined 
by resonances between the force the small halo exerts on 
the particle and the tidal force by the large halo. We as- 
sume that stripping mainly happens at primary resonance 
u(r)\ sma u = uj{R)\iarge- This leads to the following equa- 
tion for the tidal radius: 



fjXr) 
f(x R ) 



XR 



r,V m 



R S V r , 



(7) 



We take the smaller of the two estimates of r t . For 
xr > 2.2 the tidal radius is defined by the equal force 
condition. At smaller distances the orbital-internal reso- 
nance defines the tidal radius. 

Figure 1 shows tidal radius and mass within the tidal 
radius for halos at a given distance from the center of a 
group of galaxies with different masses. In this figure the 
mass of a dark halo, indicated next to each curve, is the 
virial mass of the unstripped halo, i.e. mass before the 
halo entered the cluster. As the distance from the clus- 
ter center decreases, the tidal radius of the halo and the 
mass within the tidal radius decrease accordingly. Even at 
large distances from the cluster center (R > 200/i~ 1 kpc) 
the halo radius changes significantly. For example, a halo 
with Af v i r = 10 12 /i _1 Mq and r vir = Wih-^pc at 
R = 200ft.~ 1 kpc from the center of 10 13 /i _1 M Q group 
lost only 20% of its original mass, but its radius decreased 
by a factor of two. The mass inside the tidal radius m(r t ) 
changes with R much slower m oc R°- 3 ~ a - 5 than it would 
for the isothermal distribution for which we expect m cx R. 
This is because the halos with profile Eq.(l) are more 
centrally concentrated phaio oc R~ 3 than isothermal ha- 
los (piso oc R~ 2 ). Note that the central cusp in Eq.(l) is 
not important for survival of halos at large distances: r s 
is smaller than the tidal radius r t . At smaller distances 
(R <; 2.2R S ) mass within the tidal radius decreases faster 
(m oc R) and the orbital-internal resonance defines the 
tidal stripping. It is likely that we overestimate the ef- 
fect of tidal destruction at these distances as compared to 
galaxies in real clusters because the tidal radius is small 
~ 10/i _1 kpc and the baryonic component cannot be ne- 
glected. At the same time, whether halos survive or not, 
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Fig. 1. — Tidal radii (bottom row) and masses within the tidal radius (upper row) for halos at a given distance from the center 
of a galaxy cluster of mass IO 13 /;," 1 M Q (a), 10 14 /i _1 M Q (b), and 10 15 fr _1 M Q (c). The density profiles for both the clusters 
and the halos are given by Eq.(l) with an appropriate concentration C(M). In the figures the mass of a dark halo, indicated next 
to each curve, is the mass M v j r of the halo before tidal stripping, when the halo was outside the cluster virial radius. 



they already have lost a very large fraction of their mass 
(~ 90%, exact number depends on parameters of the halo) 
when they get to R <; R s . 

2.3. Effects of energy dissipation by baryons 

Baryonic matter can loose energy by emitting radiation. 
This dissipation allows baryons to sink onto the centers of 
halos and produce a dense central region inside parent dark 
matter halo. Because a denser halo is more resistant to the 
tidal disruption, baryon dissipation clearly helps galaxies 
to survive in clusters. We consider two effects due to the 
energy dissipation, (i) The shape of the density profile 
in the central part of the halo changes without changing 



much overall structure of the halo, (ii) Global parame- 
ters of the halo (such as its characteristic radius R s and 
maximum rotational velocity V max ) change in reaction to 
the motion of baryons into the center. Another possible ef- 
fect is orbit circularization due to formation of rotationally 
supported baryonic disk. The effect of this last process, 
however, is difficult to estimate, but it likely affects only 
baryons. Below we discuss the first two effects. 

(i) The shape of the density profile in the halo center. 
The main question here is how much our estimates of the 
tidal radius would change if we assumed a steeper central 
cusp. 
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Fig. 2. — Tidal radii of halos of different masses and cir- 
cular velocities inside a cluster of virial mass 10 15 h _1 M . 
Isothermal halos (solid curves) have larger tidal radii as com- 
pared with the NFW halos (dashed curves), which means that 
isothermal halos are more stable. Both isothermal and NFW 
halos have the same maximum circular velocity. The dotted 
curve shows distance at which the tidal radius is equal to the 
characteristic radius R s - On the left from the curve the halos 
are tidally destroyed no matter what resolution is used. While 
the energy dissipation and subsequent baryonic infall increase 
the tidal radius, realistic dissipation is hardly a dramatic effect 
for halos. 



To get an estimate of the effect, we assume that sinking 
of the baryons to the halo center produces flat rotation 
curve: V c ^constant, p oc r~ 2 . For simplicity, the rota- 
tional velocity is assumed to be equal to the maximum of 
the circular velocity of the DM halo before the dissipation. 
Figure 2 shows tidal radii of halos of different masses and 
circular velocities inside a cluster of mass 10 15 ft. _1 M Q . 
Isothermal halos have larger tidal radii, which means that 
they are more stable. Nevertheless, the difference with 
the NFW halos is very small at distances from cluster 
center larger than ~ 50ft.~ 1 kpc. Even at (30 — 50)/i _1 kpc 
the the difference is only 10%-20%. The dotted curve in 
this plot shows the distance at which the tidal radius is 
equal to the characteristic radius R s of the satellite halos. 
At this distance the tidal radius of the isothermal halos 
is <J 5% larger than tidal radius of the NFW halos. Our 
TV-body simulations presented in the next section indicate 
that halos which come that close to the center of a cluster 
will be tidally destroyed regardless of what resolution is 
used. They may leave a very small leftover, which, even 
if present, is so small that it cannot represent a galaxy. 
Thus, the change in the shape of the central density profile 
cannot significantly affect survival of halos inside clusters. 

(ii) Increase of the central density due to dissipation. 
The next question is what effect for halo survival can have 



possible increase of the central density (and hence maxi- 
mum rotational velocity) due to the baryonic dissipation. 
Again, we assume that the halo after the dissipation has a 
flat rotation curve and, thus, its tidal radius can be roughly 
estimated using isothermal density profile. In Figure 2 
the halo "moves" from the NFW curve (dashed line) to a 
higher isothermal curve (solid line). For example, a DM 
halo with the virial mass of 10 u /j _:L Mpc with tidal radius 
of 7/i~ 1 kpc at lOOft^kpc from the cluster center can in- 
crease its tidal radius to 10/i~ 1 kpc, if its circular velocity 
increases from 100 km s -1 to 150 km s _1 . This may save 
the halo from the tidal destruction because its tidal radius 
is now about twice larger than its R s . Unfortunately, there 
are limits on the increase of the rotational velocity. The 
gas possesses non-zero angular momentum, characterized 
by the dimcnsionless spin parameter A (Binney & Tremain 
1987) and sooner or later becomes rotationally supported. 
For a typical value of A = 0.05, the adiabatical infall model 
(e.g., Mo, Mao, & White 1998) predicts an increase of the 
maximum rotational velocity of V c di S k/V c h a io — 1.3. (We 
use eqs. 33-34 in Mo, Mao, & White with the fraction of 
mass in disk m<j = 0.05; a correction is made to convert 
from V200 to V m ax)- Using a different approach to treat the 
baryonic infall, Avila-Reese, Firmani, & Hernandez (1998) 
arrive at the same correction: rotational velocity increases 
by a factor 1.3-1.4 for typical halo parameters. Figure 2 
shows that such increase would correspond to ~ 20 — 40% 
increase in the tidal radius. While the energy dissipation 
and subsequent baryonic infall increase the tidal radius, it 
is hardly a dramatic effect for realistic halos. 

2.4. Tidal disruption of halos: simple numerical 
simulations 

Up to this point we have treated the tidal stripping in a 
rather simplified way. Such simplified treatment is useful, 
because it gives a rough approximation for a complicated 
process. In reality (even within the framework of dark 
matter dynamics), the situation is more complicated. This 
is especially true for halos that loose a significant fraction 
of mass, in which case stability of the halos against the 
two-body relaxation is an issue. When a halo looses mass, 
the mass is lost from peripheral parts. If the trajectories 
of particles are not circular (which is typically the case), 
the central region of the halo will start to "feel" the loss on 
a dynamical time-scale. Some of particles leave the center 
and are not replaced by tidally-stripped particles; this de- 
creases the central density and leads to the expansion of 
the halo and to further loss of mass. The cycle may re- 
peat, destroying eventually the entire halo. Whether this 
process leads to halo destruction or not depends on the 
orbit of the halo and its tidal radius. 

In order to study the tidal stripping in a more realistic 
way, we run a set of small N— body simulations in which a 
halo of a few thousand self-gravitating particles moves in 
a rigid potential of a cluster. The cluster of the virial mass 
2 x 10 14 /i _1 M , a typical cluster mass for large-scale sim- 
ulations presented later in the paper, is assumed to have 
the NFW density profile eq.(l). To be consistent with our 
large simulations, we also use a ACDM cosmological model 
with the following parameters: fio = 1— = 0.3, h — 0.7, 
as = 1-0. The halo mass is chosen to be 10 12 /i _1 M Q . 
This halo has the characteristic radius of R s = 19.5/i _1 kpc 
and the maximum circular velocity of 190 kms -1 . 
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Fig. 3. — The evolution of the total bound mass (top panel) 
and the bound mass inside central 20/i _1 kpc radius (middle 
and bottom panels) for a halo with initial mass 10 12 /i _1 M 
moving inside a cluster of virial mass 2 x 10 14 /i _1 M© . Dif- 
ferent curves correspond to different force and/or mass res- 
olutions used in the experiment (each curve is marked with 
corresponding resolution) . Both the cluster and (initially) the 
halo have the NFW density profiles. The simulations were 
done with different mass mi and force e resolution. The halo 
looses 95% of its mass, but given sufficient resolution it sur- 
vives for more then lOGyr. 



Within the radius of -R200 = 163/i _1 kpc the halo has over- 
density 200 relative to the critical density of the Universe. 
At the initial moment the particles of the halo were dis- 
tributed inside -R200 m 100 equally spaced spherical shells 
in such a way that the density profile obeys eq.(l). Eq.(10) 
was used to find one dimensional velocity dispersion for 
each shell. The velocity dispersion was used to assign ve- 
locities to individual particles by throwing three random 
gaussian numbers for each velocity vector. The velocity 
distribution was isotropic. This procedure generates a sys- 
tem with the NFW profile which is almost in equilibrium. 
The shot noise (the finite number of particles) results in 
residual deviations from the equilibrium. The finite extent 
of the system also produces transient effects at the outer 
boundary. By running an isolated system with 3000 par- 
ticles for a long time (five billion years) we tested that the 
system is really close to a stationary NFW halo. We found 
only one deviation. The density distribution at the outer 
boundary (within 30% of -R200) w & s slightly smeared out. 



-300 




x(h ^pc) 



Fig. 4. — Positions of bound particles for halo presented in 
Figure 3. The time (marked next to the halo) is given in bil- 
lions of years. The solid circle indicates position of the cluster 
center. The halo initially had 952 particles and has lost most 
of its mass after 10 Gyrs. Nevertheless, even after lOGyr, 
it can still be identified as a compact centrally concentrated 
bound lump of particles. 



With the exception of the shot noise, no deviations were 
found inside 70% of i?200- 

We use the Aarseth's N-body presented in Binney & 
Trcmaine (1987). The code was modified to include the 
acceleration from the cluster. The code uses Plummer 
softening of the gravitational potential. The softening is 
thus parameterized by a softening length e. It should be 
noted that the Plummer model gives a softer force as com- 
pared with the ART code, which we used for large-scale 
simulations. We estimate that the difference in resolution 
is about factor of two: 10% error in acceleration is reached 
at ~ 2 cells in ART code as compared with ~ 3.7e in a 
code with the Plummer softening (Kravtsov et al. 1997). 

As a typical example, we have chosen an elliptical or- 
bit with the minimum distance to the cluster center of 
200ft. -1 kpc and ratio of maximum to minimum distances of 
3.35. Initially the halo was placed at distance 500/i _1 kpc 
and was given velocity 710 km s _1 (circular velocity at 
that distance) in the direction of 45 degrees to the ra- 
dius. The period of the resulting orbit is 2 x 10 9 yrs. The 
tidal radius at pericenter was 37/i _1 kpc w 1.9i? s and mass 
within the tidal radius was 1/3 of the halo's initial virial 
mass virial mass. We study the effects of the force and 
mass resolution on the evolution of the halo varying the 
resolutions by a factor of 20. 

Figure 3 shows the evolution of the total bound 
mass (top panel) and the bound mass inside cen- 
tral 20/i~ 1 kpc (middle and bottom panels). We show 
the results only if more than 15 bound particles are 
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Fig. 5. — Mass evolution for halos on different orbits. While 
halos on the highly eccentric orbits, with pericenters as small 
as 125/i _1 kpc, have survived for at least lOGyr, a halo on a 
circular orbit with radius 300fo -1 kpc was tidally destroyed af- 
ter about 7 Gyrs regardless of the used resolution. Note that 
improving mass resolution by a factor of three allows us to 
track the halo for a longer time. But results indicate conver- 
gence at mass fractions above 5%. 

found in the central region. The halo is considered "lost" 
if the number of bound particles is less than 15. For halos 
in the middle panel we fixed the mass resolution and stud- 
ied the effect of the force resolution. If the force resolution 
is not sufficient, the halo is lost. For example, with the 
force resolution of e = 10/i _1 kpc sw R s /2 the halo is lost 
after 4Gyr. Increasing the force resolution helps the halo 
to survive, but once a threshold of e = (1 — 2)/i _1 kpc is 
reached, additional increase in resolution does not change 
the mass of bound particles, which indicates convergence 
of the results. Extra resolution may even reduce the bound 
mass because of excessive two-body scattering, which was 
observed in some of our simulations. 

The bottom panel of Figure 3 illustrates the effects of 
the mass resolution. In this case the force resolution was 
fixed to e = l/i _1 kpc, sufficient for halo survival with mass 
resolution mi = 1.05 x 10 9 /i _1 Mq of the previous plot. 
The halo with three times worse mass resolution was lost 
after 6.5Gyr. As in the case of varying force resolution, 
the increase of mass resolution above some limit does not 
result in increase of the bound mass, which again points 
to convergence. This can be seen more clearly in the top 
panel: when the mass resolution is improved by a factor 
of 6 (from the dotted to the solid curve) , the mass of bound 




100 



R(kpc/h) 

Fig. 6. — The evolution of density (top) and velocity (bot- 
tom) profiles for a halo of mass 10 12 /i _1 Mq orbiting inside a 
2 x 10 14 /i _1 Mq cluster. In both panels, the solid curves show 
initial profile; the long dashed curve in the top panel shows the 
NFW profile for the initial distribution; short-dashed, dotted, 
and dot- dashed curves show profiles for halos on different or- 
bits after 10 Gyrs of orbiting in the cluster. The orbits of the 
halos are the same as the two top curves in panels of Figure 
5 (one orbit was shown twice). The curves correspond to dif- 
ferent initial to final mass ratios (shown in the legend). Note 
that the central ~ 10 kpc region is less affected by the tidal 
stripping, but when the halo looses more than 90% of its mass 
even the center becomes affected noticeably. 

particles changes only by 10%-15%. 

It should be stressed that the mass loss of the halo was 
really dramatic: after lOGyr the halo was 20 times less 
massive than at the beginning. The mass of the whole 
halo is a function of time, which steeply declines in the 
beginning and tends to level out at later moments. The 
mass loss from the central region behaves as a step func- 
tion: after every passage near the cluster center the mass 
drops by 20% - 30%. Nevertheless, the halo has survived 
(i.e., was detectable) for at least lOGyr. Figure 4 illus- 
trates the dramatic mass loss of the halo. In this plot we 
show only bound particles for a halo which initially had 
952 particles. While the halo has lost most of its mass, 
even after 10 Gyr it is still a compact centrally concen- 
trated dense lump of particles. 

Survival or destruction of halos depends on parameters 
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of their orbits in the cluster. To study the orbit depen- 
dence, we have run similar simulations varying orbit ec- 
centricities and radii. We have found that in general halos 
on circular orbits suffer larger mass loss than halos on the 
eccentric orbits. Although the pericenter of the latter may 
be smaller, the halos spend only a small fraction of the or- 
bital period in the central dense region of cluster. Figure 5 
shows mass evolution for halos on different orbits. While 
halos on the highly eccentric orbits, with pericenters as 
small as 125/i _1 kpc, have survived for at least lOGyr, a 
halo on a circular orbit with a radius of 300/i _1 kpc was 
tidally destroyed after w 7 Gyrs, regardless of the used 
resolution. 

Finally, Figure 6 shows evolution of the density and cir- 
cular velocity (V c (r) = \/GM(< r)/r) profiles of the halo 
for the same orbits as shown in Figure 5 during a 10 Gyr 
period of orbital evolution. Different orbits lead to differ- 
ent mass loss rates and the figure shows the cases in which 
the ratio of initial to final bound halo mass is 3.5, 7.7, and 
20. The figure shows that significant mass loss results in 
dramatic changes of the density and velocity profiles at 
large radii. The changes, however, are not as dramatic for 
the central halo regions. Nevertheless, in the case of ex- 
treme mass loss (« 95%) shown by the dot-dashed curves, 
both density and maximum circular velocity decrease by 
a factor of two from their initial values. 

2.5. Erasure of substructure via dynamical friction 

The dynamical friction is another effect that contributes 
to the erasure of the substructure. This is not a numerical 
effect, but by driving galaxies to the cluster center where 
they will be tidally destroyed, it can enhance numerical 
effects. The dynamical friction time for a small halo with 
mass m moving on a circular orbit of radius R around the 
large halo can be estimated using the Chandrasekhar's for- 
mula (Binneyfe Tremaine 1987) with assumptions of equi- 
librium and a Maxwellian isotropic distribution of veloc- 
ities of the DM particles. The rate of the orbital radius 
decay due to the dynamical friction is given by 
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For the density profile Eq.(l) we have d\np/d\nR = 
— (1 + 3x)/(l + x), where x — R/R s . The last of the equa- 
tions (8) can be rewritten as an equation for the velocity 
dispersion: 

do 2 



1 + 3x 2 _ GM vir 
~dx~ ~ x{l + x) <Jr ~ ~R a f(C) 
The solution of the equation is 
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Here oo is the ID velocity dispersion at R = R s . The ve- 
locity dispersion oy has a maximum oy w oo at R w 0.8i? s . 
It declines on smaller and larger radii, but the maximum 
is flat: oy w 0.78o"o &t R — 0.1i? s and oy w 0.69o- at 
R = fvir- Equations (8) and (10) define the dynamical 
friction time, which is presented in Figure 7 for differ- 
ent masses of clusters and halos. For a given halo, the 
dynamical friction time decreases as the halo moves into 
the cluster because the density of cluster increases. Note, 
however, that this decrease is countered by the halo mass 
decrease due to the tidal stripping. The dynamical fric- 
tion time is thus a varying quantity which depends on the 
cluster mass, distance from the cluster center, and halo's 
gravitationally bound mass. 

We should note that equation 8 should be considered 
as a rough estimate of the effect for at least two reasons. 
First, we assume for simplicity that orbits are circular, 
which, of course, is a very simplified view of the real halo 
orbits in clusters. This simplification, however, makes it 
possible to estimate all the interesting quantities analyti- 
cally and thus greatly facilitates computations of dynam- 
ical evolution. There is very little data as to what kinds 
of orbits are to be expected (see, however, recent work of 
van den Bosch et al. 1998). Recent numerical studies by 
Tormen (1997), Tormen et al (1998), and Ghigna et al. 
(1998) give somewhat different accounts on the distribu- 
tion of the orbital parameters. Tormen (1997) concludes 
that orbits are neither radial nor circular and that the av- 
erage eccentricity of the orbits is e « 0.5, while orbits in 
cluster analyzed by Ghigna et al. (1998) are reported to 
be mostly radial. It is not clear at present what causes the 
difference. 

Importance of the dynamical friction for a particular 
halo will depend on the halo's orbit. Halos on circular 
or low-eccentricity orbits with large radius will experience 
little gravitational drag and would never come close to the 
center and get tidally disrupted. On an eccentric orbit, 
some of these objects would have a chance to pass though 
the dense cluster core where the tidal effects discussed in 
the previous section (and gravitational drag too) would be 
considerably stronger. Note, however, that the time spent 
by a halo on highly eccentric orbit near the orbit pericen- 
ter is rather small. The strong tidal force in the cluster 
center will shock the halo during pericenter passages, but 
for most of its orbit the halo will exist in less dense (and 
thus more favorable for its survival) environments. The 
worst evolution scenario can be envisioned for a halo on 
a low-eccentricity orbit with a small radius. Such halo 
will suffer both strong and continuous tidal losses and fast 
decay of the orbital radius due to the dynamical friction. 

The validity of the Chandrasekhar's formula in spheri- 
cally symmetric, finite-size systems was also a subject of 
extensive debate over the years. One obvious drawback 
of this equation is that it docs not predict a drag expe- 
rienced by a satellite on an orbit outside of the spherical 
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system (e.g., Lin & Tremaine 1983). Most studies, how- 
ever, indicate that the equation works remarkably well. 
Lin & Tremaine (1983), for example, have addressed the 
subject numerically and have found that eq. 8 works 
well for satellites of mass <J 0.05 of the primary. The 
same conclusion was reached in other numerical studies 
(Bontekoe & van Albada 1987; Zaritsky & White 1988; 
or more recently, Cora et al. 1997). Analytical study 
by Weinberg (1986), in which Chandrasekhar's formula 
was rederived for spherically symmetric potential, also 
shows that this formalism is applicable in most cases. 
Recent analytical studies by Bekenstein & Maoz (1992), 
Maoz (1993), and Dominguez-Tenreiro & Gomez-Flechoso 
(1998) used fluctuation-dissipation approach to derive the 
energy losses due to the gravitational drag. Despite the 
more sophisticated treatment and inclusion of the various 
important details (satellite's size, not negligible mass of 
the background particles, inner velocity structure of satel- 
lite), it appears that Chandrasekhar's formula gives a good 
approximation for the energy losses. 

Dominguez-Tenreiro & Gomez-Flechoso (1998), for ex- 
ample, have taken into account both the finite size of the 
satellite and the internal velocity structure of the satellite. 
Their analysis shows that if internal velocity dispersion 
of the satellite is <; 0.5 of the velocity dispersion of the 
primary system (always the case in our study of massive 
clusters) the energy losses predicted by Chandrasekhar's 
formula are accurate to within ~ 30%. 

There is, of course, a considerably larger uncertainty 
caused by the choice of the Coulomb logarithm, because 
there is no general prescription of how to estimate the 
minimum and maximum impact parameters. The value of 
the Coulomb logarithm in clusters of galaxies is expected 
to be In A w 5 — 10. The value of In A = 8 was used in 
a recent study by Tormen et al. (1998), who compared 
dynamics of dark matter halos in clusters in the cosmo- 
logical context. They showed that even when In A is kept 
constant the Chandrasekhar's formula works remarkably 
well, predicting quite accurately decay of the satellite's 
orbital radius for the satellites as massive as <~ 0.2 — 0.5 of 
the cluster mass. This latter test is most relevant to our 
study and shows that use of the eq. 8 is justified. 

The dynamical friction time as a function of the dis- 
tance to the cluster center for a variety of cluster and 
satellite masses is shown in figure 7. In general, the results 
are hardly surprising. Dynamical friction in rich clusters 
M v ; r w 10 15 /i -1 M Q is negligible except for the most mas- 
sive galaxies near the cluster center. For poor clusters and 
groups with M v ; r <; 10 14 /i _1 M Q , the friction time is short 
as compared with the Hubble time. Adding the baryons 
would only shorten the friction time, as more mass would 
be retained within the halo's tidal radius. If a group exists 
for a sufficiently long time and does not accrete efficiently 
new satellites, the dynamical friction would produce an ob- 
ject that would look like an overmerger - a giant central 
galaxy with no other galaxies in the group 4 . The epoch 
of formation and the growth rate of groups depends on 
the parameters of a cosmological model. Thus, if resolu- 
tion is sufficient to make the numerical effects negligible, 
excessive overmerging for groups and poor clusters should 
indicate that the cosmological model is wrong. 




Distance to the center (kpc/h) 

Fig. 7. — The dynamical friction time for different masses 
of clusters and halos. For a given halo the dynamical friction 
time decreases as the halo moves into the cluster because the 
density of cluster increases. But then the halo starts to lose its 
mass, and the friction time increases again. The set of virial 
masses is the same as in Figure 1. The horizontal dashed line 
shows the Hubble time for this cosmological model. 

A matter of concern is how one could distinguish between 
numerical and physical effects as the primary cause of 
overmerging. There may be two approaches to this prob- 
lem. First, one can estimate the required resolution in the 
manner outlined in the next section and run a simulation. 
After that run another simulation with considerably higher 
resolution. Such a test may be prohibitive in terms of the 
CPU time required to run the very high-resolution run. 
However, it is probably the most certain way of testing 
for the numerical effects. High-resolution simulations of 
individual clusters are probably best suited for this kind of 
test. Indeed, a similar test was attempted by Ghigna et al. 
(1998), who analyzed two runs of the same clusters with 
different spatial resolutions. The second possible solution 
is to identify cluster's progenitor at an earlier epoch (e.g., 
z = 1) and follow evolution of a sample of test halos in the 
representative mass range. This analysis may show what 
mechanism is really at work (e.g., dynamical friction may 
be clearly observed in the orbit decay) and may allow one 
to estimate whether the numerical resolution is adequate. 



4 Notably, such systems are observed in the real Universe (Carignan et al. 1997). 



11 



2.6. What numerical resolution is required? 

Numerical experiments presented by Moore et al. (1996) 
show that halos become unstable and get quickly disrupted 
if tidal radius is smaller than ^2 — 3 times the core radius. 
The force softening, r so ft, usually has effect of creating an 
artificial core in the density distribution at scales <J r so f t . 
Thus, it can be assumed (optimistically) that halos get 
disrupted if r t <^ 2r so f t . The details of the tidal stripping 
will, of course, depend non-trivially on the mass and or- 
bital parameters of halo. However, tidal radius for a halo 
can be approximately defined (Johnston 1998; Ghigna et 
al. 1998) as the tidal radius defined by equations 6,7 at the 
pericenter of the halo orbit. Numerical studies of satellite 
dynamics (e.g., Johnston, Hcrnquist & Bolte 1996) indi- 
cate also that ~ 10 — 30% mass loss can be expected at 
every pericentric passage, which agrees well with results of 
our numerical experiments. This will have an additional 
effect if halo orbits long enough to make several pericentric 
passages. 

Thus, with an adequate mass resolution, dark matter 
halos should survive (at least during the first period of 
their orbit) inside groups of galaxies as long as their peri- 
centric r t is twice or larger than the force resolution. These 
halos may lose a large fraction of their mass, but they 
should be able to retain their identity even without the 
additional baryonic dissipation. If the required force reso- 
lution cannot be achieved, the baryons need to be included 
to alleviate the problem. 

Ultimately, survival of a particular halo will depend on 
both the mass and the force resolution. In order for two- 
body evaporation to be negligible, halos must contain J> 30 
particles. This will define the mass limit for the surviv- 
ing halos even if the force resolution is sufficiently high. A 
very optimistic estimate, thus, would be that a halo should 
consist of at least ^ 20 — 30 particles and its tidal radius 
should be larger than two resolution elements. 

As an example, we consider a rather typical simula- 
tion with particle mass of 10 10 /i _1 M Q and the resolu- 
tion of 30/i _1 kpc (e.g., Gelb & Bertschinger 1994; Ma & 
Bertschinger 1995; Tormen 1997). With this resolution 
one would naively hope to find w 10 n /i -1 M Q satellites 
- the virial radius of a halo with this mass is more than 
twice larger than the force resolution. The top panel in 
Figure 1(b) indicates, however, that halos of this ini- 
tial virial mass cannot be found inside a cluster of mass 
10 14 /i _1 M Q : due to the insufficient mass resolution they 
lose so much of their mass that they are destroyed by the 
tidal force. The mass resolution of these simulations was 
sufficient to find what is left of a halo with 3 x 10 n /i _1 M 
at distances ;> 80/i _1 kpc. However, the bottom panel 
in Fig. 1(b) shows that tidal radius for galaxy-size ha- 
los is smaller than two resolution elements within cen- 
tral ~ 300/i _1 kpc and thus no such halos should be ex- 
pected there. Therefore, lack of the force resolution re- 
sulted in erasure ( "overmerging" ) of even fairly massive 
(~ lO 12 ^ 1 Mq ) halos within the SOOh^kpc. 

But even a considerably better resolution may not be 
sufficient if we deal with a really massive cluster. For ex- 
ample, Carlberg (1994) simulated a 2.2 x lO 15 /^ 1 M 
cluster with mass resolution of 2.27 x 10 9 /i _1 M Q and 
the Plummer force softening of e — 9.7/i~ 1 kpc. Accord- 
ing to Gelb & Bertschinger (1994), the effective resolution 



for the Plummer force is 2.6e. This gives the resolution 
w 25/i~ 1 kpc, which we use as a limit on the tidal radius of 
resolved halos. Analysis of tidal radii for a cluster of this 
mass shows that due to the insufficient force resolution, no 
halos should exist at distances smaller than <J 290/i~ 1 kpc. 
Halos of mass M vir < 10 n /i _1 M should be tidally de- 
stroyed at R — 590/i~ 1 kpc. This is consistent with what 
Carlberg (1994) found in his simulation. 

What resolution is required for halos to survive? The 
above examples show that the answer depends on the mass 
of the cluster and on the mass of the halo one would like 
to resolve. The force resolution should be (significantly) 
smaller than the minimal tidal radius of a halo (defined by 
the orbit's pericenter). For a 10 14 ft- _1 M cluster the tidal 
radius for a massive halo of mass Mh ;> lO 11 /^ 1 M at 
a distance of R £ (60 - 70)/i~ 1 kpc is r t w 10/i _1 kpc. The 
force resolution should probably be better than 3/i _1 kpc 
for such halo to survive in the cluster. Because the halo at 
this distance loses 80%-90% of its original virial mass, and 
because realistically one needs at least 20-30 particles to 
identify a halo, the particle mass should be smaller than 
i=a 10 9 h~ 1 M Q . These estimates may be optimistic for 
halos that orbit in cluster for several periods due to con- 
tinuing mass loss (see §2.2) and possible effects of halo-halo 
heating (Moore et al. 1996). 

To some extent, the answer is clear. In order to allow 
a galaxy-size halo to survive, the force resolution must be 
much smaller than its tidal radius and the halo must be 
represented by many particles. For a typical large galaxy 
of the mass lO 11 ^ -1 M and tidal radius of 10/i _1 kpc in 
a <~ 10 14 — lO 15 ^^ 1 M cluster, the resolution must be 
of (0.5 — 3)/i _1 kpc and particle mass - <J 10 9 /i _1 M . 
There is a number of further caveats related to analysis 
of the halo remnants. The most important of them is the 
question of whether it is possible to recover properties of 
the original halo from the properties of the remnant. Un- 
fortunately, the answer to this question is not clear and 
various parameters may be affected in different and com- 
plicated ways. Most of the the halos orbit in cluster for 
periods less then 10 Gyrs, and although their peripheral 
parts are being stripped, the evolution of the central re- 
gions is not dramatic (see Fig. 6). 

3. SIMULATIONS 

We simulate the evolution of 128 3 cold dark matter 
particles in two cosmological models: a flat low-matter 
density CDM model with cosmological constant (ACDM; 
Slo = 1 — Cl\ = 0.3; h = 0.7; cr% = 1-0) and a model with 
a mixture of cold and hot dark matter (CHDM; Qq = 1; 
Q, v = 0.2; h = 0.5; cr 8 = 0.7 ). The CHDM simulation 
followed trajectories of additional 2 x 128 3 hot particles. 
The mass fraction of hot matter, Q v , is equally split be- 
tween two types of neutrino (Primack et al. 1995). Both 
models were normalized to be consistent with COBE DMR 
observations (Bunn & White 1997). Normalization of the 
ACDM model is also consistent with observed abundance 
of galaxy clusters (Viana & Liddle 1996), while normal- 
ization of the CHDM model may be slightly higher than 
suggested by the data (Gross 1997). Both simulations are 
initialized with the same set of initial random numbers in 
order to reduce effects of the cosmic variance. 

To achieve sufficiently high mass resolution, the size of 
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the simulation box for both the ACDM and the CHDM 
models is chosen to be rather small - 15/i _1 Mpc. To 
test the effects of the box size, an additional simulation 
of 30/i~ 1 Mpc box has been run for the ACDM model. 
The mass of a cold particle in the 15/i _1 Mpc box simu- 
lations (for 30ft- _1 Mpc box particle mass is 8 times larger) 
is mi = 1.33 x lO 8 /^ 1 M for the ACDM model and is 
« 2.7 times larger for the CHDM model. If we assume that 
;> 30 particles are needed to identify a halo, we expect to 
be able to identify halos as small as 4 x 10 9 /i _1 M Q . This 
corresponds to 4 x 10 10 /i _1 M Q halo before tidal stripping, 
if 90% of mass was stripped by the cluster. 

The simulations were done using the Adaptive Refine- 
ment Tree (ART) code (Kravtsov et al. 1997). The 
code used a 256 3 uniform grid on the lowest level of res- 
olution and seven levels of refinement for the 15ft. -1 Mpc 
box. Each refinement level doubles the resolution. The 
seventh refinement level corresponds to the dynamical 
range of 32,000 and the resolution of w 0.5/i _1 kpc. The 
30/i Mpc run had six levels of refinement and its reso- 
lution is w 2/i~ 1 kpc. The code refines an individual cell 
on a given level L if the number of particles in this cell 
(as estimated by the cloud-in-cell method) exceeds some 
threshold N t h(L). The threshold is N t h — 5 for high levels 
and N t h — 10 is set for low levels L = 0, 1 (Kravtsov et 
al. 1997). This choice of thresholds ensures that refine- 
ments are introduced only in the regions of high-particle 
density and prevents the two-body relaxation effects. The 
increase in spatial resolution corresponding to each suc- 
cessive refinement level is accompanied by decrease of the 
integration time-step by a factor of 2. The simulations 
were started at z% = 30 when the rms of the density fluc- 
tuations in the simulation box was 6 « 0.27 — 0.32. 

The dynamic range of the simulation is justified by the 
following two considerations. First, the code integrates the 
evolution in comoving coordinates. Therefore, to prevent 
degradation of force resolution in physical coordinates, the 
dynamic range between the start and the end (z = 0) of 
the simulation should increase by (l + z,): i.e., for our sim- 
ulations 256 x (1 + Zi) = 7680. Second, the code reaches its 
peak resolution in the highest density regions inside virial 
radius of the DM halos. To resolve a halo, we need at least 
<~ 10 resolution elements per halo size, which justifies the 
dynamic range of ;> 10, 000. 

Particle trajectories were integrated with the step in ex- 
pansion factor of Aao = 0.0015 on the zero level uniform 
grid, and with time step Aa^ = Aa /2 L on a refinement 
level L. This gives an effective number of steps of 82,000 
on the seventh level of refinement. In physical units, the 
smallest time step at z = corresponds to 1.15 x 10 5 years. 

4. HALO IDENTIFICATION 

Finding halos in dense environments is a challenge. The 
most widely used halo-finding algorithms: the friends-of- 
friends (hereafter FOF, e.g., Davis et al. 1985) and the 
spherical overdensity algorithm (e.g., Lacey & Cole 1994; 
Klypin 1996) - are not acceptable (Gelb & Bertschinger 
1994; Summers et al. 1995). The friends-of-friends (FOF) 
algorithm merges together apparently distinct halos if link- 
ing radius is too large or misses some of halos if the radius 
is too small. Adaptive FOF (van Kampcn 1995) seems to 
work better. However, our experiments show that in prac- 



tice it is difficult to find an optimal scaling of the linking 
radius with the density for a general case. 

We have developed a version of the FOF algorithm, 
which we call "hierarchical friends-of-friends" . This algo- 
rithm uses a fixed set of hierarchical linking radii and thus 
does not have problems adaptive FOF algorithm has. The 
algorithms, cither adaptive or hierarchical, cannot work 
using only geometrical means to identify a halo. In the 
very dense environments they pick up many fake halos. 
This can be improved by taking into account dynamical 
information to decide whether a halo is real or not. The 
DENMAX algorithm (Bertschinger & Gelb 1991; Gelb & 
Bertschinger 1994) or its offspring SKID (Governato et al. 
1997) make a significant progress - they remove unbound 
particles, which is important for halos in groups and clus- 
ters. Another approach to deal with "flukes" is to check 
whether the halos in question were distinct halos at earlier 
epochs. 

Recently, Summers et al. (1995) tried to perfect the 
idea of Couchman & Carlberg (1992) to trace the history 
of halo merging and to use it for halo identification. Start- 
ing at an early epoch, Summers et al. identify halos using 
the FOF algorithm with linking radius corresponding to 
the "virial overdensity" of 200 and then trace particles 
belonging to halos at later times. It appears that it is im- 
possible to make a working algorithm. Halos interact too 
violently. A large fraction of mass is tidally stripped from 
some halos and a large fraction of mass is being accreted 
by others. However, the idea of using a set of epochs is 
too good to be abandoned. To avoid the problems with 
the "direct approach" (from past to the future), we have 
decided to try a reverse logic. Instead of asking the ques- 
tion "where is now the halo that collapsed at some earlier 
epoch" , we ask "did the halo that we find at present exist 
at an earlier time?" . Thus, we supplement our hierarchi- 
cal FOF algorithm with an algorithm which checks if halos 
existed at previous moments. 

The algorithm which finds halos as maxima of mass in- 
side spheres of a given overdensity works better than the 
plain FOF, but no fixed overdensity limit can find halos in 
both low and high density environment. However, the at- 
tractive simplicity of this algorithm makes it worth in look- 
ing for ways of improving it. We will discuss our improve- 
ments to such an algorithm in §4.2. Besides the problem 
of finding whether a halo is real, a problem of "missing" 
halos should be kept in mind. As was discussed in the 
previous section, while some halos may survive in dense 
environments, others may be destroyed due to numerical 
effects. These halos will be missing at z — and "present 
to past" approach does not help. 

Some of the problems that any halo finding algorithm 
faces are not numerical. They exist in the real Universe. 
We select a few of the most typical difficult situations. 

1. A large galaxy with a small satellite. Examples: the 
LMC and the Milky Way or the M51 system. Assuming 
that the satellite is bound, do we have to include the mass 
of the satellite in the mass of the large galaxy? If we do, 
then we count mass of the satellite twice: once when we 
find the satellite and then when we find the large galaxy. 
This docs not seem reasonable. If we do not include the 
satellite, then the mass of the large galaxy is underesti- 
mated. For example, the binding energy of a particle at 
the distance of the satellite will be wrong. The problem 
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arises when we try to assign particles to different halos in 
the effort to find masses of halos. This is very difficult 
to do for particles moving between halos. Even if a par- 
ticle at some moment has negative energy relative to one 
of the halos, it is not guaranteed that it belongs to the 
halo. The gravitational potential changes with time, and 
the particle may end up falling onto another halo. This is 
not just a precaution. This actually was found very often 
in real halos when we compared contents of halos at dif- 
ferent redshifts. Interacting halos exchange mass and lose 
mass. We try to avoid the situation: instead of assign- 
ing mass to halos, we find the maximum circular velocity, 
y/GM/ R\ m ax, which is a more meaningful quantity than 
mass from observational point of view. 

2. A satellite of a large galaxy. The previous situation 
is now viewed from a different angle. How can we estimate 
the mass or the rotational velocity of the satellite? The 
formal virial radius of the satellite is large: it coincides 
with the virial radius of the host halo. In order to find the 
outer radius of the satellite, we analyze the density profile. 
At small distances from the center of the satellite the den- 
sity steeply declines but then it flattens out and may even 
increase. This means that we reached the outer boundary 
of the satellite. We use the radius at which the density 
starts to flatten out as the first approximation for the ra- 
dius of the halo. This approximation can be improved by 
removing unbound particles and checking the steepness of 
the density profile in the outer part. 

3. Tidal stripping. This is not a numerical effect and is 
not due to a "lack of physics" . Very likely this is what hap- 
pens to real galaxies in clusters. Their peripheral parts, re- 
sponsible for extended flat rotation curves outside of clus- 
ters, are lost when the galaxies fall into a cluster. Thus, 
if an algorithm finds that 90% of mass of a halo identified 
at early epoch is lost, it does not mean that the halo was 
destroyed. This is a normal situation. What is left, given 
that it still has a large enough mass and radius, is a galaxy 
halo. 

4.1. Hierarchical friends- of -friends algorithm 

In order to find substructures at vastly different over- 
densities, we use hierarchical friends-of-friends algorithm 
(HFOF). This algorithm simply applies the FOF algorithm 
with a set of different (hierarchical) linking lengths. In our 
analysis we use 4 hierarchical levels, in which case the set 
consists of linking lengths starting from the small value 
l v ir/8 and larger values obtained by doubling this value: 

/ — f*i)X7* 1 4j I — ^viv / ^ ) cUld I ^VZT* ' AA^G Cclll \ ^VIT 

the lowest (with respect to the corresponding overdensity 
threshold) level and I — l V i r /8 the highest level. The link- 
ing length l v i r corresponds to the virial overdensity of an 
object. We assume the virial overdensity of 200 and 340 
(l mr = 0.21 and l mr = 0.170 for SCDM/CHDM and for 
ACDM models, respectively (e.g., Lahav et al. 1991; Eke 

— —1/3 

et al. 1996); here, I = n ' is the mean interparticle sepa- 
ration and no is the mean particle density in the simulation 
box. For our halos the smallest linking length, I = l v i r /8, 
corresponds to overdensity of s=s 10 5 . At each level of the 
hierarchy identified clusters of particles (halo candidates) 
are marked if none of their particles belongs to a marked 
higher-level cluster. Finding halo candidates first at the 
highest possible level is important because some of higher- 



level clusters merge into larger halos at lower levels. In 
practice, we find that in high-resolution simulations there 
is a wealth of small clumps on all levels in large cluster-size 
halos, defined on the lowest level of the hierarchy l vir . 

Another important feature of this algorithm is use of a 
particle distribution at an earlier epoch. The FOF algo- 
rithm works on a snapshot of the particle distribution and 
generally tends to identify particle clusters that are linked 
at this moment just by chance. Typical examples would 
be either a "bridge" connecting two particle clusters or 
small satellite clusters on highly eccentric orbit that move 
temporarily beyond the virial radius of a massive halo. 
Therefore, one must check the stability of every identified 
halo candidate. The easiest way to make such a check is 
to find whether a given halo candidate exists at an earlier 
moment. We perform such check by running the HFOF 
algorithm with the same set of linking lengths using posi- 
tions of particles at z = 1 and check both the existence and 
one-to-one correspondence of the progenitor particle clus- 
ters. We consider a candidate halo to be "stable" if it has 
one (or two) progenitor (s) and it is the only descendant of 
the progenitor(s). 

Details of this analysis are as follows. We select the two 
most massive progenitors of each halo and check whether 
they combined contain more than a threshold number of 
particles of the halo at z = 0. We search for these pro- 
genitors on the next lower level of hierarchy to take into 
account the fact that the size of an unevolving object at 
z = 1 doubles in comoving coordinates. The level is not 
reduced if it reaches the virial overdensity. We sum the 
two most massive progenitors in order to allow one ma- 
jor merging. In fact, there are very few cases in which 
the masses of three progenitors were of the same order of 
magnitude. We used the thresholds 70%, 50% and 30% of 
mass of the halo at z = 0. The threshold 70% is too high: 
in many cases halos accrete more than 30% of their mass 
between z = 1 and z = 0. The algorithm would fail to find 
many halos with this threshold. On the other hand, there 
is little difference in the number of identified halos using 
thresholds of 30% and 50%. However, the algorithm would 
not find the most massive and the most dense halos in the 
center of the large groups because the mass of those halos 
increases substantially between z = 1 and z = due to 
merging with small halos and accretion of single particles. 
To avoid this, we also include halos which contain more 
than a minimum number (100-500) of particles at z = 1. 
We found that the result almost does not depend on the 
chosen number because at z = 1 these progenitors contain 
already considerably more particles. 

The code checks also whether the particles found in the 
progenitor represent a substantial fraction of its mass. The 
importance of this criterion is clear from the following ex- 
ample. Close to the very massive and very dense halos 
the FOF algorithm finds many small clumps. At an ear- 
lier moment, all of these small clumps belong to the same 
progenitor of the massive cluster around which they were 
found at z = 0. Each lump, however, represents only a tiny 
fraction of the progenitor. To avoid misidentifications the 
algorithm accepts only halos whose particles represent a 
substantial fraction of the mass of that progenitor. We 
used values of 10% and 30% as thresholds for the mass 
fraction of the progenitor . Results are not sensitive to 
the particular choice of the threshold. Both criteria (total 
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Fig. 8. — A halo identification by the HFOF. Top row: Circles and the triangle represent positions of identified halos at the level 
lvir/8 (overdensity of » 10 5 ). Points without circles show fake halos (flukes). Depths of the projections are 200ft _1 kpc for the left 
panel and 5/i~ Mpc for the right panel. Bottom row: left panel shows position of dark matter particles of the fake halos at z = 1, 
while the "real halos" are shown in the bottom right panel. Most of particles of the five "real halos" are in very compact halos 
whereas the fake halos are very puffy (spread over large groups seeing on the top panel of 



mass and mass fraction of the progenitor) result in selec- 
tion of stable halos. 

An example of the halo identification in very dense en- 
vironment is presented in Figure 8. In the top left panel it 
shows the inner region of a group-size halo with the trian- 
gle denoting the position of the identified extended massive 
central "galaxy" of the group. The points surrounded by 
circles represent dark matter particles in five identified ha- 
los at the level l V i r /8 (overdensity s=s 10 5 ). Points without 
circles show fake halos (flukes) : the halos are found by the 
FOF algorithm at this level, but they do not satisfy the 
fraction of progenitor criterion. Bottom left panel shows 
position of dark matter particles of the fake halos at z = 1, 
while the "real halos" are shown in the bottom right panel. 
Most of particles of the five "real halos" are in very com- 
pact halos, whereas the particles of the fake halos are in 



extended and "puffy" configurations. 

4.2. Bound density maxima algorithm 

In addition to the hierarchical FOF, we have developed 
another halo-finding algorithm, which uses ideas of the 
DENMAX algorithm (Bertschinger & Gelb 1991; Gelb & 
Bcrtschinger 1994). Just as the DENMAX, our algorithm 
first finds positions of the density maxima on some scale 
and then removes unbound particles inside the halo radius 
(hence, the name: Bound Density Maximuma (BDM)). 
However, the algorithm finds maxima and removes un- 
bound particles in a different from the DENMAX way. 
The algorithm can work by itself or in conjunction with 
the hierarchical FOF. In the latter case, it takes positions 
of halos from the HFOF, and then removes unbound par- 
ticles and finds parameters of halos. The version of the 




Fig. 9. — An example of a poor cluster in a ACDM sim- 
ulation. At z — the cluster has a virial mass of 2 x 
10 13 /i _1 M , virial radius of 500h~ 1 kpc, and velocity dis- 
persion of ~ 500fcm/s. The figure shows all DM particles 
(~ 250, 000) in a sphere of radius 1.5/i -1 Mpc with the centered 
at the cluster. The particles are color-coded on a grey-scale 
according to the log 2 of the local density (the density is esti- 
mated as a number of particles at r < 8h~ kpc from a given 
particle). Figure 10 shows the halos identified in this volume 
by the bound density maxima algorithm. 



Fig. 10. — Dark matter halos identified by the bound density 
maxima algorithm in the volume shown in Figure 9 (volume and 
projection are the same). 



BDM code used here is available for use by astrophysical 
community (see Klypin & Holtzman 1997). 

In order to find positions of halos we choose a smoothing 
radius r sp of a sphere for which we find maxima of mass. 
This defines the scale of objects we are looking for, but 
not exact radii or masses of halos. Radius of a halo can be 
either larger or smaller than r sp . For example, if we are 
interested in galaxy-size halos, it is reasonable to choose 
r sp ~ (10 — 15) kpc. If we search for galaxy groups, an 
appropriate choice is r sp ~ (200 — 300)kpc. Then we place 
a large number of the spheres in the simulation box. The 
number of the spheres is typically an order of magnitude 
or more larger than the number of expected halos. For 
each sphere we find its center of mass and the mass inside 
r sp . The center of the sphere is displaced to the new cen- 
ter of mass and this process is iterated until convergence. 
Depending on specific parameters of the simulations, the 
number of iterations ranges from 10 to 100. This process 
finds local maxima of mass within r sp . Some of the max- 
ima will be found many times. We remove duplicates and 
keep only one halo for each maximum. Halos with too 



small number of particles (typically 5-10) and halos with 
too low central overdensity are removed from the final list. 

Figures 9 and 10 illustrate typical output of the algo- 
rithm. Figure 9 shows a group-size halo identified in our 
30/i _1 Mpc ACDM simulation. The particles, shown in the 
figure inside the sphere of radius 1.5/i _1 Mpc centered on 
the group, are color-coded on a grey-scale according to the 
logarithm of the local density, estimated as a number of 
particles within 8/i _1 kpc of the particle. Figure 10 shows 
dark matter halos identified by the BDM algorithm in this 
volume. All of the halos, that can be seen in Figure 9 as 
tight dark clumps of particles, are identified by the halo 
finder. 

Once centers of potential halos are found, we start the 
procedure of removing unbound particles and finding the 
structure of halos. We place concentric spherical shells 
around each center. For each shell we find mass of the dark 
matter particles, mean velocity, and the velocity disper- 
sion relative to the mean. In order to determine whether 
a particle is bound or not, we estimate the escape veloc- 
ity (Eqs.5) at the position of the particle. If velocity of a 
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Fig. 11. — Examples of profiles of small halos in the ACDM model. Each column of plots corresponds to the same halo. Mass 
and radius of each halo (within virial overdensity 340) are shown in the top panels. The dashed curves are for the halos before the 
removal of unbound particles, and the full curves are for bound particles only. The top row of panels shows the circular velocity 
V c i rcu iar = (G M (R) / R) 1 ^ . The middle row shows velocity dispersion of the dark matter particles, and the bottom row presents 
the overdensity profiles. The right halo is an example of an isolated halo in which most particles are bound to the halo. The two 
halos on the left show examples of small satellite close to a large halo. In the central ~ 20/i _1 kpc part of both halos the velocity 
dispersions are small and almost constant. But at 30/i~ 1 kpc the velocity dispersion starts a dramatic increase indicating presence 
of a massive object. 



particles is larger than the escape velocity, it is assumed 
to be unbound. We estimate the maximum rotational ve- 
locity V max and radius of the maximum r max = 2r s using 
the density profile for the halo. Because V max and r max 
must be found before the unbound particles are removed 
and because the mean velocity is also found using all par- 
ticles (bound and unbound), the whole procedure can not 
be done in one step. We start by artificially increasing the 
value of the escape velocity by a factor of three. Only par- 
ticles above the limit are removed. We find new density 
profile, new mean velocities, new V max and r max . The es- 



cape velocity is again increased, but this time by a smaller 
factor. The procedure is repeated 6 times. The last itera- 
tion does not have any extra factors for the escape velocity. 

Removal of unbound particles is crucial in the case when 
a halo with a small internal velocity dispersion moves in- 
side a large group. For example, if a halo with a circular 
velocity of 100 km s _1 moves with velocity 500 km s _1 in- 
side a group, a dark matter particle bound to the group has 
kinetic energy 25 times larger than the kinetic energy of a 
particle of the halo. Even if only l/10th of particles within 
the halo radius belong to the group, the whole halo will 
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Fig. 12. — The same as in Figure 11 but for medium and large halos. In this case the masses are given for constant outer radius 
of 122ft" J kpc. The mean overdensity at that radius for all halos is significantly larger than the virial value. But the maximum of 
rotational velocity is well within the distance. 



have positive energy and will be treated as fake. Removal 
of unbound particles salvages the halo even if the real halo 
particles constitute as little as 1/4 of the total number of 
particles within the halo radius. This estimate is valid in 
the case of a compact halo moving through homogeneous 
field of high velocity particles of the group. The situa- 
tion is worse if the particles of the high velocity field are 
very lumpy. The worst case is the collision of two equal 
mass halos. At the moment when the halos overlap, the 
code does not find any bound component - both halos are 
missed. The chance of such event is very small because the 
distance between centers of halos should be smaller than 
- (10- 15)/i- 1 kpc. 

The effect of the unbound particle removal on the halo 



profiles is illustrated in Figures 11 and 12. The halos were 
identified in the simulation of the 30ft _1 Mpc ACDM sim- 
ulation. In the Figure 11 the right column shows circular 
velocity, velocity dispersion, and density profiles of a "nor- 
mal" halo with a small fraction of unbound particles, pri- 
marily in the outer regions. The middle and left columns 
show profiles of small satellite halos located inside or close 
to a massive halo. In the central ~ 20ft~ 1 kpc of both 
halos the velocity dispersions are small and almost con- 
stant. The circular velocities are about what one should 
expect for these values of the velocity dispersions. But at 
30ft _1 kpc the velocity dispersion increases dramatically, 
indicating a presence of a dense background of fast moving 
particles belonging to the massive halo. Figure 12 shows 
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typical examples of medium- and high-mass halos. 



4.3. Comparison of the HFOF and BDM algorithms 

The goal of both of the described algorithms is to find 
positions of stable halos in a given simulation. Neverthe- 
less, the two algorithms are quite different. The HFOF, 
for example, computes mass, spin parameter, angular mo- 
mentum, shape and total binding energy of the halos with- 
out removal of unbound particles. The algorithm, how- 
ever, does not assume spherical symmetry in these calcu- 
lations. The BDM algorithm computes properties of the 
halos (e.g., mass, velocity, density profile) after removing 
unbound particles, the procedure based on the assumption 
of spherical symmetry. Therefore, one must expect that 
the algorithms may compute slightly different properties 
for the same halos. 

Nevertheless, the most important information deter- 
mined by halo finders is positions of DM halos. Our tests 
show that the algorithms find exactly the same halos, if 
the halos contain more than a couple hundred particles. 
The agreement is about 95 % for halos with more than 
50 particles and about 90 % for halos with more than 30 
particles. The small differences in the threshold criteria 
for the selecting halos account for the 10 % differences at 
the low mass end of the halo distribution. But overall, 
the agreement is very good. For example, the five small 
and the big central halo in a high density region shown 
on the top left panel of Figure 8 have been found by both 
algorithms. We also compare halos in a statistical way. In 
Figure 17 we compare the correlation functions of halos 
found by the two algorithms. For this comparison we have 
selected halos with rather low limit in maximum circular 
velocity V max > 90 km/s. In case of the ACDM simula- 
tion (lower panel) the two correlation functions coincide 
within the expected scatter due to the statistical noise. 
This indicates that the thresholds are chosen in an equiv- 
alent manner. In case of the CHDM simulation (upper 
panel) the correlation functions differ more. Very likely 
this happened because of the slight mismatch in the mass 
limits. As expected, the two correlation functions differ 
mainly on scales less than 100 kpc due to different thresh- 
olds used in the algorithms. 

5. RESULTS 

5.1. Luminosity function of galaxies and M/L in groups 

The observed tight correlation between the 21-cm line 
width W ~ 1V C i rc and infrared luminosities for spiral 
galaxies (e.g., Aaronson & Mould 1983; Bureau, Mould, 
& Staveley-Smith 1996; Willick et al. 1996; Giovanelli 
et al. 1997) can be used to estimate luminosities of 
galaxy-size halos in A-body simulations (an alternative 
method for assigning luminosities to the DM halos and 
constructing the luminosity function was proposed re- 
cently by Roukema et al. 1997). This is probably the best 
one can do when dealing with dissipationless simulations. 
Unfortunately, the observed Tully-Fisher relation is not 
defined as accurately as one would hope for. This is es- 
pecially true for low-luminosity galaxies. In the following, 
we apply the empirical Tully-Fisher relation determined 
for galaxies with absolute blue magnitudes rag < —15. 
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Fig. 13. — Comparison of the luminosity function of galaxy- 
size halos in the LCDM simulations with the CfA data (dashed 
curve). The solid circles show results for the 15ft _1 Mpc box. 
The results for 30ft -1 Mpc box are shown with open circles. 
The first bin of the 30ft" 1 Mpc simulation at M = —21.45 has 
8 galaxies; the bin at M = —17.8 has 350 galaxies. There are 
146 galaxies brighter than Bj = —20 and 1340 brighter than 
Bj = — 18 in the 30ft~ Mpc box simulation. 

It should be kept in mind, however, that there may be sig- 
nificant systematic deviations from the relation used here 
for galaxies with magnitudes of raj > —17. 

Elliptical galaxies pose another problem. The Faber- 
Jackson relation indicates that velocity dispersion (and 
thus the dark matter mass) can be used to estimate galaxy 
luminosity. The relation, however, is not very tight. In 
principle, by tracing the merging history of each halo, 
we can make more realistic estimates of the star forma- 
tion rates and the luminosities. Here, we will use a sim- 
ple but reasonable prescription; we assume that an el- 
liptical galaxy is ~ 1 magnitude dimmer than a spiral 
galaxy with the same maximum circular velocity. This 
assumption is motivated by the fact that mass-to-light 
ratio of elliptical is 2.5-3 times higher ratio than that 
of sprials spirals. It is likely that the fraction of ellip- 
ticals is significant only at the the high and low mass 
ends of mass function. We assume therefore that all 
halos with V C i rc > 350 km s" 1 and half of the halos 
with V C i rc < 80 km s _1 host elliptical galaxies. Halos 
in the "grey area" V c i r c = (200 — 350) km s _1 have a 
gradually increasing probability to host an an elliptical: 
cx 1 / 3(V/ 200km /s) 2 . 

It is important to stress the use of maximum circu- 
lar velocity and/or velocity dispersion of halo as a mass 
indicator. Problems of finding dark matter halos and 
determining their masses in numerical simulations high- 
lighted in §4, make it virtually impossible to make a 
meaningful assignment of mass to a halo in a group 
or a cluster. On the other hand, velocity dispersion 
or maximum circular velocity can be determined more 
or less reliably with even moderate particle statistics. 
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Fig. 14.— The same as in Figure 13, but for the CHDM 
15/i -1 Mpc box simulation. 

We believe therefore that use of the maximum circular ve- 
locity as a proxy for halo mass is more attractive than use 
of potentially unreliable and biased mass estimates. 

In order to estimate luminosities of "galaxies" hosted 
by DM halos using their maximum circular velocities, 
we use the following Tully-Fisher relation in the I-band: 
Mi - Uog(h) = -21.0 - 6.8(logW - 2.5). The slope of 
the relation is as given from Willick et al. (1996) for field 
spirals, while the zero point was adopted from Giovanclli 
et al. (1997). The I-band magnitudes were shifted to the 
blue magnitudes as M B = M/ + 1.5 (Pierce & Tully, 1992). 
Because circular velocities in observations are never mea- 
sured at large galactocentric distances, we set an upper 
limit of r max < 50/i _1 kpc for the radius of the maximum 
of the rotational velocity. 

The main caveat of the above luminosity assignment 
scheme is that we use the maximum circular velocities of 
halos, discarding the disk contribution. However, as we 
have discussed in §2.3, this may result in maximum er- 
ror of 30% in the circular velocity. Also, Figure 6 shows 
that maximum circular velocity of halos orbiting in clus- 
ters may decrease by ~ 20 — 50% due to the tidal stripping. 
The luminosity assignment may be somewhat biased for 
such halos if baryons are not stripped as efficiently as the 
dark matter. 

Figure 13 shows the luminosity function (LF) of galaxy- 
size halos identified in the two of our ACDM simula- 
tions. The solid circles show results for the 15ft. _1 Mpc 
box. There is a significant tail of low luminosity galaxies 
(Bj > —16.5), which matches the corresponding tail of 
the CfA luminosity function (Marzke et al. 1994) and the 
luminosity function of cluster galaxies (Smith et al. 1997; 
Lopez-Cruz et al. 1997). The simulation of the larger 
box of 30/i~ 1 Mpc (open circles) has insufficient mass res- 
olution to probe this low-mass tail, but the two LFs are 
consistent in the region of overlap. We believe that the 
LFs in Figure 13 are most reliable in the magnitude range 
of Bj = —17 — 20. In this range the number of of ha- 



los in each luminosity bin is sufficiently high to make the 
poisson errors insignificant and results do not depend on 
the assumed fraction of ellipticals or on the maximum al- 
lowed radius for the circular velocity. The first bin of the 
30/i _1 Mpc simulation at M = —21.45 contains 8 halos, 
while the bin at M = —17.8 contains 350 halos. The lu- 
minosity function in the ACDM model in this range of 
magnitudes is systematically lower by a factor of 1.5-2 
than the CfA luminosity function. This is consistent with 
deeper samples, which give lower normalization for the lu- 
minosity function (e.g., Loveday et al. 1992). 

The luminosity function of galaxy-size halos in the 
CHDM simulation is significantly higher than both the 
ACDM and the CfA luminosity functions. With the same 
set of parameters as for the ACDM model, the LF in the 
CHDM model is a factor of 4-5 higher than the luminos- 
ity function in the CfA catalog. In order to reconcile the 
model with the observational data, we reduced the limit on 
the radius for the rotational velocity to r max = 25/i _1 kpc 
and raised the fraction and the limiting magnitude for 
small elliptical galaxies. Figure 14 shows the CHDM lu- 
minosity function that best matches the CfA luminosity 
function. It is still systematically higher than the CfA 
LF, but it might be acceptable due to small volume of the 
simulation. 



5.2. Halo dynamics in groups: velocity bias, M/L, and 
constrains on £1 

The observed mass-to-light ratios of galaxy groups and 
clusters, (M/L) b » 150-400 (e.g., Bahcall, Lubin & Dor- 
man 1995), are often used as an argument in favor of the 
low-density universe with VLq w 0.2 — 0.3. We have used 
the halos identified in poor galaxy clusters and groups in 
our simulations to study their dynamical properties and 
estimate the mass-to-light ratio of these clusters using the 
same prescription to assign luminosity to halos as was used 
in the previous section. Figures 15 and 16 present dif- 
ferent properties of groups of galaxies in the simulations. 
Centers of the groups were found using search radius of 
r S p — O-250/i^ 1 Mpc and no removal of unbound particles 
was done in this case. Radius of groups was estimated 
at the overdensity limit of 200 for the CHDM model and 
of 340 for the ACDM model. The results clearly indicate 
that the M/L ratio increases with the mass of the group. 
However, there is an indication that in the ACDM simu- 
lation M/L ratio for massive groups flattens at the level 
of w 3OO/i _1 (M /i ). Note, that models with n = 0.3 
(ACDM) and f2 = 1 (CHDM) reproduce the observed 
M/L ratios equally well, although due to the small vol- 
ume, we can only probe masses of <; 3 x 10 13 ft. _1 M Q in 
the CHDM simulation. It appears thus that mass-to-light 
ratio of galaxy groups of mass <; 3 x 10 13 /i -1 M Q is not 
a very good indicator of the total matter density in the 
universe. 

Several authors have suggested a possible existence 
of the velocity bias b v = a halo /o~ dm'- systematic differ- 
ence between rms velocities of galaxies and dark mat- 
ter particles (Carlberg, Couchman & Thomas 1990; 
Carlberg 1994; Colafrancesco, Antonuccio-Delogu & 
Del Popolo 1995). The existence of the velocity 
bias would have impact on the determination of clus- 
ter masses using galaxy dynamics and other analyses. 
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Fig. 15. — Properties of groups in the ACDM simulations as 
a function of the total mass of the group. V rms ,dm is the root- 
mean-square velocity of dark matter particles in the group. 
The ratio of the rms velocities of galaxies in the group to 
V rms ^dm is shown in the second panel. Results for groups with 
more than 3 satellites are shown. In small groups with mass 
less than 10 13 ft _1 Mq there is on average a significant ve- 
locity bias, which is likely related to the dynamical friction. 
There is no indication of the velocity bias for larger groups. 
The number of galaxies with V C irc > 50 km s _1 in a group is 
shown in the third panel from the bottom. The mass-to-light 
ratio of groups (in units of h^ 1 Mq /Lq^b) is shown in the 
top panel. 

Figures 15 and 16 show that there is no evidence for strong 
velocity bias for halos in groups in our simulations; on av- 
erage, for objects ;> 1O 13 /i -1 M0 the rms velocities of 
galaxies and dark matter particles are equal. Although 
theoretical arguments suggest that a mild velocity bias 
should exist (due, for example, to effects of dynamical 
friction), the uncertainties of current simulations favor 
absence of the velocity bias (b v = 1) and excludes values 
of the bias b v <; 0.8, allowing, however, for a possible 
existence of a mild bias b v w 0.8 — 0.9. Note that for small 
groups of galaxies (lO 12 ft _1 M £ M mr £ lO 13 h _1 M ) 
a significant velocity bias is observed. This bias could 
probably be attributed to the strong dynamical friction 
effects operating in these systems. 

5.3. Small scale two-point halo correlation function 

Analyses of recently completed galaxy surveys resulted 
in a very accurate determination of the galaxy two-point 
correlation function, £(r), at the scales of ;> 20 — 50/i _1 kpc 
(e.g., Baugh 1996). The comparison of the observed £(r) 
with the two-point correlation function of mass in dissipa- 
tionless A-body simulations has shown that an antibias of 
galaxies is required at small scales in order to reconcile the 
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Fig. 16. — Properties of groups in the CHDM simulation. 

models with observations (Klypin, Primack & Holtzman 
1996; Jenkins et al. 1998). In order to check whether 
any significant antibias exists for the dark matter halos 
(which could be associated with observed galaxies), we 
have constructed the halo-halo two-point correlation func- 
tion for the halos found in our simulations using both 
hierarchical friend-of-friend and bound density maxima 
algorithms. Figures 17 and 18 show correlation functions 
of halos and dark matter in the simulations. The corre- 
lation function is clearly affected at scales ^> 500ft.~ 1 kpc 
by the finite box size (as shown by comparison of £(r) 
for small and large box simulations of the ACDM model) . 
However, it is interesting to examine the relative behavior 
of the dark matter and halo correlation functions at scales 
<; 500/i _1 kpc. At very small scales (less than 100/i _1 kpc) 
the galaxies are more clustered (biased) relative to the 
dark matter; at larger scales the effect is the opposite: 
galaxies are slightly antibiased. This is observed in all 
simulations and it is valid for all limits on masses of 
galaxies. It is interesting that the antibias of the mag- 
nitude 0.7-0.9 seen in the Figure 17 is almost exactly 
what is needed for the ACDM model to be compatible 
with observational data on the power spectrum in the 
range of wavenumbers k = (0.1 — l)ftMpc -1 (Klypin et 
al. 1996; Smith et al. 1998; Jenkins et al. 1998). The 
antibias of halos at these small scales is very likely re- 
lated to the tidal destruction and dynamical friction of 
halos in groups. The 100/i _1 kpc — l/i _1 Mpc range is the 
range where we expect both of these processes to work. 
This conjecture seems to be confirmed by results of larger 
simulations (Colin et al. 1998; Kravtsov & Klypin 1999). 
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Fig. 17. — Correlation functions of the dark matter (long 
dashed curves) and galaxy-size halos in the ACDM (bottom) 
and the CHDM simulations (top) in the 15/i~ Mpc boxes. 
Galaxies with rotational velocity larger than 90 km s were 
selected. Short-dashed curves are for galaxies identified with 
the hierarchical friends-of-friends algorithm. The solid curves 
are for halos identified using bound density maxima algorithm. 
The circles are for halos identified using the bound density 
maxima algorithm in the 30/i _1 Mpc box simulation with the 
same limit on the rotational velocity. 



6. CONCLUSIONS 

We have presented arguments that the overmerging 
problem is mostly due to inability of a numerical code to 
provide a sufficient numerical resolution to prevent tidal 
destruction of galaxy-size halos by the tidal forces of group 
or cluster and have estimated the resolution needed to pre- 
vent such disruption. We argue that although energy dissi- 
pation by the baryonic component helps galaxies to survive 
in clusters, at distances J> (50 — 70)/i _1 kpc from the clus- 
ter center the gravity of the dark matter alone is enough 
to keep them alive. 

The main result of this work is estimate of the nu- 
merical resolution needed to overcome the overmerging 
problem. The results of our analytic estimates and nu- 
merical experiments show that although it is feasible to 
overcome overmerging in pure A^-body simulations, resolu- 
tion required to avoid artificial destruction of galaxy-size 
halos (mass ;> lO n /i _1 M ) is quite high. For viable 
CDM models and realistic halo profiles this resolution is 
<; 2/i -1 kpc in force and <; 10 9 /i _1 M© in mass. This re- 
quires simulations of > 10 particles with dynamic range 
of 10 in spatial resolution for statistically significant 
cosmological volumes (<~ lOOMpc), which remains chal- 
lenging with the current computers and numerical codes. 
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Fig. 18. — Dependence of the galaxy correlation function 
on mass in the ACDM simulation of 30/i -1 Mpc box. The cor- 
relation function increases with the rotational velocity, but all 
curves show the same tendency: positive bias on small scales, 
slight antibias on (100 — 1000)/i _1 kpc scales, and no bias on 
larger scales. Absolute values of the correlation functions are 
affected by the finite box size. 

This makes simulations focused on the individual clusters 
of the type presented in Ghigna et al. (1998) a viable 
alternative. 

Unfortunately, even in the case of sufficient resolution, 
when halos do survive, the identification of the DM ha- 
los in the cores of galaxy groups and clusters in purely 
dissipationlcss simulations remains a challenge. In the en- 
vironments that dense, most of halo's dark matter will 
be be tidally stripped, which makes it difficult to iden- 
tify the leftover on the very dense, smooth background of 
high-velocity dark matter particles streaming around and 
through the halo. We have presented two new halo find- 
ing algorithms designed to identify satellite halos located 
inside the virial radius of a more massive host halo: the hi- 
erarchical friends-of-friends and bound density maxima al- 
gorithms. Both of our algorithms find practically the same 
halos, which are stable (existed at previous moments) and 
gravitationally bound. 

To exploit the fact that overmerging is (at least to a cer- 
tain degree) "beaten" in our simulation, we consider sev- 
eral statistics of galaxy-size halos in our simulations and 
compare them to the corresponding observed statistics of 
galaxies. We use a simple scheme, based on the empirical 
Tully-Fisher relation, to assign a luminosity to the DM ha- 
los. The luminosity function of "galaxies" (i.e., galaxy-size 
halos assigned a luminosity) in the ACDM model repro- 
duces the luminosity function of the CfA catalog (Marzkc 
et al. 1994) reasonably well. Both the simulations and the 
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CfA catalog have an upturn in the number of faint galax- 
ies (m B > —17). However, magnitudes of faint "galaxies" 
in the simulations rely on a highly uncertain extrapolation 
of the Tully-Fisher relation and on uncertain assumption 
about the fraction of elliptical galaxies at these magni- 
tudes. The number of "galaxies" predicted by the CHDM 
simulation is significantly higher than in the case of the 
ACDM simulation with the same initial random numbers. 
We failed to produce as nice a fit to the observational data 
as for the ACDM simulation. At this stage it is difficult to 
judge if this is a significant problem for the model or not. 
Due to the small size of our simulation boxes, one may 
argue that simulations with a large box will tend to pro- 
duce lower luminosity function keeping at the same time 
the M/L of galaxy groups intact. Larger simulations are 
needed to clarify the situation. 

The mass-to-light ratios of galaxy groups in the simula- 
tions ~ (200 — A00)h^ 1 match those observed reasonably 
well. It was argued (e.g., Bahcall, Lubin, & Dorman 1995) 
that dynamics of galaxy groups favors the low-Jl Universe. 
Our results show that mass-to-light ratios of groups of 
mass < 3 x 10 13 /i -1 M Q is insensitive to the matter den- 
sity. The halos in the CHDM model are clustered more 
strongly than the dark matter and one cannot save the 
argument for a low-f2 Universe by assuming that groups 
in the CHDM model have too large fraction of galaxies. It 
seems that groups occupy too small fraction of the volume 
and thus their M/L ratios are not representative for the 
Universe as a whole. 

Comparison of the halo and matter correlation functions 
indicates that halos are antibiased on 100 kpc - 1 Mpc 
scales. The antibias of the magnitude 0.7-0.9 found in the 
simulations is needed for the ACDM model to be compat- 
ible with observational data on the power spectrum in the 
range of wavenumbers k — (0.1 — l)/iMpc _1 (Klypin et al. 
1996; Smith et al. 1998). We attribute the antibias to the 
dynamical friction in groups of galaxies. The friction tends 
to drag some galaxies to the very central part of groups 
where they merge the central galaxy and disappear (see 



Kravtsov & Klypin 1999 for a more detailed analysis). 

Results of this paper can be used in design of the fu- 
ture numerical simulations. We have shown that efficient 
halo finding algorithms can be developed to identify grav- 
itationally bound satellite halos inside the virial radius of 
the other halos. Our analytical estimates and numerical 
experiments show that the numerical resolution required 
to overcome the overmerging, although quite high, can be 
achieved with current numerical codes and computer hard- 
ware. The main challenge is thus purely computational. 
This is also true for the simulations that include dissipative 
hydrodynamics; while alleviating or obliterating some of 
the problems of dissipationless simulations, they are com- 
putationally more intensive. Both numerical approaches 
have a number of caveats and potential biases, which could 
only be avoided with inclusion of more realistic physics. 
The latter appears to be unavoidable, because we can- 
not reliably predict observed galaxy properties (and hence 
mimick the selection criteria of the observational catalogs) 
without realistic physics. Fast increase in computational 
capability of modern computers and recent developments 
of new efficient numerical algorithms make the perspective 
for advances in this direction look good. 

We thank Joel Primack for careful reading of the 
manuscript and comments, Simon White for discussions, 
and anonymous referee for constructive criticism and 
useful suggestions which have helped to improve con- 
tent and presentation of the paper. We are grateful to 
Avishai Dekel for providing us with computer resources 
at the Hebrew University. This work was funded by the 
NSF and NASA grants to NMSU, and the collaborative 
NATO grant CRG 972148. SG acknowledges support 
from Deutsche Akademie der Naturforscher Leopoldina 
with means of the Bundesministerium fur Bildung und 
Forschung grant LPD 1996. Our simulations were done 
at the National Center for Supercomputing Applications 
(Urbana-Champaign, Illinois) and on a Power Challenge 
supercomputer at Hebrew University. 



REFERENCES 



Aaronson, M., & Mould, J. 1983, ApJ, 265, 1 

Avila-Reese, V., Firmani, C, & Hernandez, X., 1998, ApJ, accepted 

Bahcall, N., Lubin, L., Dorman, V. 1995, ApJ, 447, L81 

Baugh, CM. 1996, MNRAS 280, 267 

Bekenstein, J. D., Maoz, E. 1992, ApJ, 390, 79 

Bertschinger, E., & Gelb, J.M. 1991, Comp. Phys., 5, 164 

Binney, J., & Tremaine, S. 1987, "Galactic Dynamics" (Princeton: 

Princeton University Press) 
Bontekoe, Tj. R., van Albada, T. S. 1987, MNRAS, 224,, 349 
Brighenti, F., & Mathews, W.G. 1997, tistro-ph/9707070 
Bunn, E.F., & White, M., 1997, ApJ, 48U7Ti 

Bureau, M., Mould, J.R., & Staveley-Smith, L. 1996, ApJ, 463, 60 
Carignan, C, & Freeman, K. C. 1988, ApJ, 332, L33 
Carignan, C, Cote, S., Freeman, K.C., &; Quinn, P.J. 1997, AJ, 114, 
1313 

Carlberg, R.G, Couchman, H.M.P., Thomas, P.A. 1990, ApJ 352, 
L29 

Carlberg, R.G. 1994, ApJ, 433, 468 

Carrol, S.M., Press, W.H., & Turner, E.L. 1992, ARA&A, 30, 499 
Colin, P., Klypin, A. A. , Kravtsov, A V j.Khokhlov, A.M. 1998, ApJ 



submitted (preprint lastro-ph/9809202 ) 
Cora, S. A., Muzzio, j. C, Vergne, M. M. 1997, MNRAS, 189, 253 
Colafrancesco, S., Antonuccio-Delogu, V., Del Popolo, A. 1995, ApJ 

455, 32 

Couchman, H.P.M., & Carlberg, R ApT W) AW 

Courteau, St., & Rix, H.-W. 1997, 
Davis, M., Efstathiou, G., Frenk, C 
292, 371 



istro-ph/9707290 
S., fc White, S.D 



M. 1985, ApJ, 



de Blok, W.J.G., McGaugh, S.S. 1997, MNRAS, 290, 533 
Dominguez-Tenreiro, R., & Gomez-Flechoso, M.A. 1998, MNRAS, 
in press 

Eke, V.R., Cole, S., & Frenk, C.S. 1996, MNRAS, 282, 263 

Evrard, A.E. 1997, MNRAS, 292, 289 

Faber, S.M., & Gallagher, J.S. 1979, ARA&A, 17, 135 

Forman, C, Jones, C, & Tucker, W. 1985, ApJ, 293, 102 

Gelb, J., & Bertschinger, E. 1994, ApJ, 436, 467 

Ghigna, S., Moore, B., Governato, B., Lake, G., Quinn, Th., Stadel 

J. 1998, MNRAS, 300, 146 
Giovanclli, R., Haynes, M., da Costa, L., Freudling, W., Salzer, J., 

& Wegner, G. 1997, ApJ, 477, LI 
Governato, F., Moore, B., Cen, R., Stadel, J., Lake, G., & Quinn, 

Th. 1997, New Astronomy, 2, 91 
Gross, M.A.K., 1997, PhD Thpsis, T.TC Santa Cni7 (available 



electronically at tittp : //f ozzie .gsf c .nasa. gov/index . html ) 
Jenkins, A. et al., iyy8, ApJ, 4y9, 2U 
Johnston K.V., 1998, ApJ, 495, 297 

Johnston.K.V., Hcrnquist, L., & Bolte, M. 1996, ApJ, 465, 278 
Katz, N., Hernquist, L., Weinberg, D.H., 1992, ApJ 399, L109 
Klypin, A.A., Primack, J., Holtzman, J. 1996, ApJ, 466, 13 
Klypin, A. A. 1996, in "Dark Matter in the Universe", p. 419, eds. 
S. Bonometto, J. Primack, A. Provenzale, IOS Press, Amsterdam, 
Oxford, Tokyo, Washington DC 
Klypin, 



A. A., & Holtzman, I 1QQ7 prnprint ag -t-T-A- r h/Q71?017 (thcr-q 



de 



is available at ittp : //astro . nmsu . edu/ aKiypin/ pmcode . html ) 
Kravtsov, A.V., 6i Klypm, A. A. lyyy, m preparation 



23 



Kravtsov, A.V., Klypin, A. A., & Khokhlov, A.M. 1997, Astrophys. 

J. Suppl., Ill, 73 
Lacey, C, Cole, S. 1994, MNRAS, 271, 676 

Lahav, O., Lilje, P.B., Primack, J.R., & Rees, M.J. 1991, MNRAS, 
251, 128L 

Lin, D.N.C., & Tremaine, S. 1983, ApJ, 264, 364 

Lopez-Cruz, O., Yee, H.K.C., Brown, J. P., Jones, C, & Forman, W. 

1997, ApJ, 475, 97L 
Loveday, J., Peterson, B.A., Efstathiou, G., Maddox, S.J., 1992, 

ApJ, 390, 338 
Ma, C.-P. & Bertschinger, E. 1995, ApJ, 434, L5 
Maoz, E. 1993, MNRAS, 263, 75 

Martimbeau, N., Carignan, C, & Roy, J.-R. 1994, AJ, 107, 543 
Marzke, R., Geller, M., Huchra, J., & Corwin, H. 1994, AJ, 108, 437 
Merritt, D. 1985, ApJ, 289, 18 

Mo, H.J., Mao, S., White, S.D.M., 1998, MNRAS, 295, 319 
Moore, B., Katz, N., & Lake, G. 1996, ApJ, 457, 455 
Navarro, J., Frenk, C, & White, S.D.M. 1995, MNRAS, 275, 720 
Navarro, J., Frenk, C, & White, S.D.M. 1996, ApJ, 462, 563 (NFW) 
Navarro, J., Frenk, C, & White, S.D.M. 1997, ApJ, 490, 493 
Persic, M., Salucci, M., & Stel, F. 1996, MNRAS, 281, 27 
Pierce, M., & Tully, B. 1992, ApJ, 387, 47 

Primack, J., Holtzman, J., Klypin, A., & Caldwell, D. 1995, Phys. 

Rev. Lett., 74, 2160 
Rix, H, W, 1997, Proceedings of the DM 1996 workshop, Sesto, 

astro-ph/9611040 



Juinn, P.J., Rocca-Volmerange, B. 

Thonnard, N. 1985, ApJ, 
Holtzman, J. 1998, 



Roukema, B.F., Peterson, B.A. 

1997, MNRAS 292, 835 
Rubin, V.C., Burstein, D., Ford, W.K.. 

289, 81 

Smith, C, Klypin, A., Gross, M., Primack, J. 

MNRAS, 297, 910 
Smith, R., Driver, S., U Phillipps, S. 1997, MNRAS, 287, 415 
Summers, F.J., Davis, M., & Evrard, A. 1995, ApJ, 454, 1 
Tormen, G. 1997, MNRAS, 290, 411 

T T m. n a. n; g w; n A., Syer, D., 1998, MNRAS submitted 



_astro-ph/9712222| ) 

cteE BosctT F ' . Leads. 



G E , L 



ke, G., Stadel, J. 1998, ApJ, 



in press (preprint astro-ph/981122£ ) 
van Kampen, E. 1995, MNRAS, 273, l 295 
Viana, P.T.P., & Liddle, A. 1996, MNRAS, 281, 323 
White, S.D.M. 1976, MNRAS, 177, 717 
Weinberg, M.D. 1994a, AJ, 108, 1398 
Weinberg, M.D. 1994b, AJ, 108, 1403 
Weinberg, M.D. 1997, ApJ, 478, 435 

White, S.D.M., Navarro, J.F., Evrard, A.E., Frenk, C.S. 1993, 
Nature, 366, 429 

Willick, J., Courteau, St., Faber, S., Burstein, D., Dekel, A., k. 

Kolatt, T.S. 1996, ApJ, 457, 460 
Zaritsky, D., & White, S.D.M. 1994, ApJ, 435, 599 



